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Abstract 

We give an exponentially-accurate normal form for a Lagrangian particle moving in a 
rotating shallow- water system in the semi-geostrophic limit, which describes the motion in the 
region of an exponentially-accurate slow manifold (a region of phase space for which dynamics 
on the fast scale are exponentially small in the Rossby number). We show how this result 
is related to the variational asymptotics approach of |Oli05| : the difference being that on the 
Hamiltonian side it is possible to obtain strong bounds on the growth of fast motion away 
from (but near to) the slow manifold. 

Our normal form approach extends to numerical approximations via backward error anal- 
ysis, and extends to particle methods for the shallow-water equations, where the result shows 
that particles stay close to balance over long times in the semi-geostrophic limit. 

1 Introduction 

Semi-geostrophic approximations are obtained from the rotating shallow- water equations in the low 
Rossby and Burger number limit. In |Oli05| . a systematic approach for deriving these models is set 
out; based on looking for a near-identity change of coordinates for the solution domain which makes 
the Lagrangian affine (linear in velocity), resulting in a equation which only has slow dynamics. At 
each order in the Rossby number, there is some choice of parameters available, and this choice can 
be crucial in obtaining a well-posed PDE. After the change of coordinates, the equations of motion 
can be obtained by applying Hamilton's variational principle to the approximated Lagrangian. 

In this paper we use exponentially-accurate normal form theory to investigate the rotating 
shallow- water equations in the limit of low Rossby and low Burger number. Exponentially-accurate 
normal form theory was first popularised through the work of |Nek77| (very clearly reviewed 
and reformulated in |Pos99p on perturbed intcgrable Hamiltonian ODEs. Also important in the 
literature, |Nei81j gave an exponentially-accurate estimate for systems with rapidly rotating phase, 
and |B(l(T87j developed estimates for highly oscillatory mechanical systems. A first application of 
exponentially-accurate normal form transformations in the context of geophysical fluid dynamics 
can be found in |WS00j . The principle analytic obstruction to a result in our case is the requirement 
that the free surface height remains sufficiently smooth; hence, we restrict our result to individual 
particle trajectories moving in a fluid with a specifled free surface height (a problem which has 
applications in tracking advected tracer particles: see the summary for details), to toy ODE systems 
and to numerical solutions obtained with the Hamiltonian particle-mesh method for the rotating 
shallow-water equations. The central tool of the normal form theory is the choice of a suitable 
symplectic transformation for the particle positions and momenta (this transformation may be 
more general than the cotangent-lifted coordinate transformations used in |Oli05j ). resulting in a 
Hamiltonian system where the small coupling between fast and slow time-scales is more explicit. 
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More explicitly, if the model equations are solved with initial conditions which are close to balance, 
then they will stay there for very long time intervals. We also note that our results are more general 
than the type provided in |Wir04| , where the existence of a slow manifold up to exponentially small 
terms is shown. 

The rest of the paper is organised as follows. In section [51 we introduce the model system of a 
particle moving in a rapidly-rotating shallow-water flow, and discuss the ideas behind geostrophy 
and semi-geostrophic approximations. In section 12.11 we introduce the variational structure of 
these equations from both of the Lagrangian and Hamiltonian viewpoints, and in section [2.21 we 
discuss higher-order balanced models in the context of the variational structure. This background 
then allows us to introduce the exponentially-accurate normal form theory in section 12.31 and 
we show how this theory extends to numerical solutions obtained with symplectic time-stepping 
methods in section 1531 Section applies this theory to particle methods for the rotating shallow- 
water equations. Section [3.11 discusses geostrophic balance for these methods in the context of 
exponentially-accurate normal form theory, whilst section 13.21 illustrates the theory with some 
numerical results. Finally there is a summary and outlook in sectional 



2 A model system 

The shallow-water equations (SWEs) on an /-plane are 

= +.fv - gfix, (1) 
= -/-" - gf^y, (2) 

= -/i (Ua; + Vy) . (3) 



Du 
Dv 



Here fj, = fj,(t,x,y) is the fluid depth, g is the gravitational constant. / is twice the (constant) 
angular velocity of the reference plane. 



D 
Dt 



{.) = {.\ + u{\^ + vn, (4) 



is the Lagrangian or material time derivative, and subscripts denote partial differentiation with 
respect to that variable. See |Sal99j for a derivation of the SWEs and their relevance to geophysical 
fluid dynamics. 

In non-dimensionaliscd variables, the SWEs can be written in the form 

Du B , , 

^ Dv B 



Dfi 

In 



-pL {Ux + Vy) , (7) 



where Ro = U/fL is the Rossby number and B = [Ln/L)^ is the Burger number. Furthermore, 
Lr = y/gH/ f is the Rossby radius of deformation, [/ is a typical advection velocity, L is a typical 
length scale, and H is the mean fluid depth for the problem under consideration. Semi-geostrophic 
theory is concerned with the limit of (0-0 with Ro = 0{e), B = 0(e) as e — > 0. For simplicity, 
we use Ro = B = e from now on. See jSal99j for a derivation of non-dimensionalised equations and 
the semi-geostrophic scaling limit. 

Instead of investigating the full SWEs, we start in this paper with the following simpler model 
problem. We assume that the layer depth ^{et, x, y) is a given function of time t and space x = 
(x, ■ We may then investigate the motion of a single fluid parcel with coordinates q ~ [q^, Qy)^ ■ 
The corresponding Newtonian equations of motion arc given by 

d'^qx , dqy , ^ X /DA 

= - thA£t,qx,qy), (8) 



= -M9„(e^,fe,9y), (9) 

We next rescale time and introduce the new time-scale t = et. We denote time derivatives with 
respect to r by overdot, e.g., dq^/dr = (jx- We also write the second order equations JH))-® as a 
system of first order equations by introducing the momentum p = {px,Py)'^ , i-S-, 

Px = +Py - efig^{T,qx,qy), (10) 

Pv = -Px - £ f^qy{T,qx,qy), (11) 
Qx = Px, (12) 

qy = Py. (13) 

These equations may be expressed in a more compact form; 

P = ^2P - e Vq/x(r,q), J2 = { q) , (14) 



(15) 



which is the formulation wc will work with in this paper. A discussion of the relationship between 
the Hamiltonian structure of this system and the Hamiltonian structure for the full SWEs is given 

in unnn. 

The intuitive idea behind the semi-geostrophic approximation |Sal99| is that the solutions of 
C3J"CSJ consist of inertial oscillations with period Tj = 2it = ©(e") and slow geostrophically 
balanced motion on a (slow) time-scale Tq = C>{e). Furthermore, intuition suggests that the 
inertial oscillations arc primarily governed by the linear equations 

P = 'hp, (16) 
q = P, (17) 

while the slow, geostrophically balanced, parcel dynamics is characterised by the reduced (nonlin- 
ear) system 

= J2p-£VqAi(T,q), (18) 
q = p, (19) 

or, equivalcntly, 

q = -eJ2VqAi(T,q). (20) 
2.1 Hamiltonian and variational formulations 

We start with a Hamiltonian formulation of the system H14() - (|15|l . As a further simplification we 
assume that the layer-depth fi is time-independent and introduce the potential energy function 
y(q) := /-t(q) to make the link to classical mechanics more transparent. 

We rewrite the equations (|14|I -H15 |I using a non-canonical symplectic structure operator 



so that 

with Hamiltonian 



z = JV,iJo(z), (22) 
1 



Ho{z)^K{p) + eV{ci), /i(p) = -p^p, (23) 

and phase space variable z = (p"'",q^)^ G R^. 

Another approach is to rewrite the system of first-order equations as a second-order equation 



q - Jaq + £VqF(q) = 0, 



(24) 



and to note that H24|l is the Euler-Lagrange equation for the Lagrangian functional 



dr 



(25) 



The reduced model (|20|) with = V is also Hamiltonian with phase space q G M^, structure 
matrix jj, and Hamiltonian function Hg{q) = eV{q). The associated Lagrangian functional is 
given by 



Note that ^ differs from ^ by the missing kinetic energy term ||qjp/2. 



(26) 



2.2 Higher order balance and the 'semi-geostrophic' approximation 

Clearly, one would like to derive more sophisticated reduced models; a systematic approach has 
recently been given in jOli05j . In this section we summarize the main results from |Oli05| before 
we develop our novel approach in section [2.31 

All known derivations of balanced (semi-geostrophic) models start from the assumption that 



q = -£J2Vqy(q) + 0{e^) = 0{e). 



(27) 



Under this assumption, we may formally collect terms of equal order in e and rewrite (|25() as 

C = s J dT[Lo + eLi], where Lo ^ - V (q) , Li = ^||qf . (28) 

We now introduce a coordinate transformation i/'s : ^ M'^ and transformed coordinates q^ via 



q = V',(q,) = q, + eF,{q,) + 0{e^). 
Following |Oli05| (but note our different notation), we set 



Fi(qe) = -^Vqy(qe) + AVqV^(q,) = -^J2qe + AVqF(qe) + 0(e). 



Hence we obtain a Lagrangian function 

Ci=e j dT [Zo + eLi+0(£2)] 
in the transformed variable q^, where 



1 



V^V{q,fj2qe-\\\V^V{q,)\\ 



Upon dropping 0{e^) terms in %M\. we finally obtain the transformed Lagrangian 
1 



dT 



-j^qe J2qe - eV{q,)+e 



VqU(qe)^ J2qe-e'A||VqU(qe)|| 



(29) 



(30) 



(31) 



(32) 



(33) 



In the context of this paper, the choice A = —1/2 is of particular interest. Note that this 
choice corresponds to Salmon's large-scale semi-geostrophic approximation (see |Sal99l lOliOSp . 
The associated reduced equations of motion can be derived from the Lagrangian functional 



C 



Isg 



dT 



^||Vqy(qe)|p + iqJj2q.-eV^(qe) 



The associated Euler-Lagrange equations are explicitly given by 



J2qe - eVqU(qe) + y Vq||VqU(qe)|l2 = 0. 



(34) 



(35) 



These equations are canonical with Hamiltonian 

i7i.g(qe) = eV{<i,) - ^||Vqy(qe)f , (36) 

and structure matrix Jigg = jj. 

A number of questions arises at this stage. 

1. The assumption (|27(l is essential for the derivation of the reduced Lagrangian H^lfl . However, 
even if (|27|l holds at the initial time, how can one ensure that it holds true for solutions of 
the full equations (|24|l over finite (or infinite) time intervals? 

2. Following |Oli05j . one can formally generalize the coordinate transformation to obtain higher- 
order balanced equations. Can this process be continued to arbitrary precision and what 
rigorous error bounds in terms of e can be obtained? 

3. What role is played by the parameter A? 

We will provide an answer to these questions from the perspective of Hamiltonian perturbation 
theory in the following sections. We will first restrict the discussion to the simple model system 
(|14|I - H15|1 before considering a generalization to particle methods for the SWEs in section |3| 

We also wish to point the reader to the reports |GO05al IGO05b| , which attempt to answer the 
above questions based on higher-order coordinate transformations 

q, = q + eFi + e^F2 + ■■■+ e'"F„, + 0(£™+i) (37) 

and associated higher-order expansions 

C^^e J dT[Lo + eLi + ---+ e"Z^ + 0{e"'+^)] (38) 

of the Lagrangian functional for m > 1. 



2.3 Exponentially accurate normal forms 

In this section we view the finite dimensional system ()22ll from the Hamiltonian side and look for 
symplectic coordinate transformations to a normal form. 

The precise aim of the normal form transformation is to find a symplectic (with respect to the 
structure matrix J) change of coordinates, which transforms the system to a form in which the fast 
and slow variables are completely separated up to a small remainder term. As the transformation 
is symplectic, the equations of motion can be obtained simply by substituting the change of coor- 
dinates into the Hamiltonian and writing down Hamilton's equations. This normal form strategy 
is well-developed and we only summarise the main steps. See, e.g., jNei84| for the first derivation 
of an exponential estimate and |Cot04j for an application to systems of type ()22|l. 

Our aim is to find a near-identity change of coordinates so that 

Hn = Hoo^„ = K + eGn + e"+'i?„ , (39) 

where 

{G„,if} = 0, (40) 

with {•, •} being the Poisson bracket obtained from J. 
We define 5'„+i recursively by writing 

*n+l = *n O i , *0 = H, (41) 

where $i.en+iF„+i is the time-1 flow map of the Hamiltonian vector field produced by the Hamil- 
tonian hmction H = £"+^F„+i. Recursive substitution gives 

= Hn o $i,e"+if„+i - A- + eGn + e"+\Rn + {K, K+i}) + 0(£"+'), (42) 



and so we need to choose Fn+i so that 

Rn + {K,F„+i}^Rr., (43) 

where 

{Rn,K}^0. (44) 

Equation (|43|l is known as the homological equation. The solution -Fn+i, and the non-resonant part 
Rn of Rn can be calculated via 

1 

Rn ^ 7^ drRnO^r^K, (45) 

1 '■^ 



Fn+1 = 7f dTT{Rn-Rn)o<^r,K, (46) 

where ^r,/^ is the time-r flow-map of the Hamiltonian vector field produced by K , which is time- 
periodic with period T = 2tt. To close the recursion we define 

G„+i = G„ + £"i?„ (47) 

and terms of size ©(e""*'^) and smaller are collected in a new residual function Rn+i- 

When the potential V is analytic and bounded on an appropriate compact set in the complex 
plane, we can choose to make the number of iterations scale with 1/e and obtain an exponentially- 
accurate normal form as summarised in the following theorem: 

Theorem 2.1. Let he some compact subset of phase space containing the solution trajectory, 
and let V be analytic in BrJ(^ ( the union of complex balls of radius r with centres in J^). Then 
there exists a near-identity change of coordinates '^^ such that 

H = Hoo^J^ = K + eG + R, (48) 

where 

{G,A-} = 0, (49) 
G and R have e-independent bounds on J(f and c > is some constant. 

Proof. Essentially the method is to obtain a Cauchy estimate for each remainder and to 
truncate the above algorithm after n steps where n scales with 1/e. Defining the norm on functions 

sup \f{z)\, (50) 
z£B,^,jfr 

we obtain iterative estimates for Rn and Fn+i from equations (|45I46(I |P5s99j : 

I-Rnir < |i?„|r, \Fn+i\r < 27T\R„\r ■ (51) 

Simple adding and subtracting of terms gives 

iJ„o$i_,„+ij.„^^ = A-+[eG„ + e"+MA-,F„+i}+e"+ii?„] 
+K o $i,,.+i;^„^, ~K~ e"+i{A', F„+i} 

+ e{Gn O *l,e" + ii^„ + i - G„) 

+ e"+l(i?„ o $i,,„ + i;^„^, - Rn). (52) 

The term in the square brackets is eG„+i, and it remains to estimate this term, and the remainder. 
If we choose S so that i — {n + 1)6 > and i7„ is analytic in B^-nS^ then 

\Gn+l\r-{n+l)5 = |G„ + e"i?„ |r-(n+l)(5i 

< |G„|r-Ti(5 + e"|-Rri|r-Ti(5) (53) 



having used the estimate for i?„. Estimates for the remainder Rn+i can be obtained using the 
mean- value theorem combined with a Cauchy estimate for the gradient. To estimate the second 
hne in equation l|52() note that 

= e"+i / {/C,i^„+i}o<D,.,„+ij.,,^^dt-e"+i{/C,F„+i}, 
Jo 



Jo 

[ (Rn-Rn)o'i't,er^ + lF„+,-{Rn-Rn)dt, (54) 
JO 

which may be bounded using the mean-value theorem and a Cauchy estimate by 

|/C O $1 ,„ + ij,^^^^ - /C - e"+M^>.9n+l}|r-(n+l)5 

< e"+i / |(i?„-i?„)°$t,."+ij.„+, -(i?„-i?„)|,„(„+i),dt, 







< e^^"+^^|V(i?„ — Rn)\r-{n+l}5\Fn+l\r-(n+l)S 

^ J \Rn ^ Rn\r-ns\Rn\r~{n+l)57 

27re2("+l) 

< ^ \Rn\r-n5\Rn\r-nS- (55) 

Similar estimates give 

"+iF„ + i G',i|r_(ri+l)5 < „ IGrilr-jiil-Rrilr— n5: (56) 

1 27re2("+i) 

e" l-Rn O $l,e" + i_F„+i — -Rn|r-(ri-Hl)(5 < j |^ri |r-n(S | |r-n(S • (57) 

Suppose we have: 

n 

|G„|r-„5 < Co — j , |i?„|r-«5 < Co (^yj , (58) 



where cq is a positive constant, and ci 47rco, then we get 

n.+l 

V S J 

1=1 

from equation H53|l . and 



|G„+i|r-(„+i)5 < Co ^ ' ^^^^ 



27re"+i 27r 

l-^n+l |r— (n+l)(5 ^ | -^n |r— n(5 | -^n | r — n(5 ^~ | l^-— 7^(5 [t- — 7^(5 

27re"+l 

H ^ |-^n |r — n5 \Rn \r — n5 ; 

, (I)" (1^.. I (-) '), (.0, 

from equations H55I57|I . If 



then |G„+i|r_(„i)5 < co, and 



These estimates are true for all n by induction if we set cq = \V\r. Finally wc choose n to scale 
with e from the formula 

" 2r " 



(63) 

ecie_ 

where the square brackets indicate rounding up to the nearest integer, and set 5 = r /2n. Equation 
(|61|l can then be satisfied, and we get 

e"+'|i?n|r/2 < eco(^)" < eco ^ ' = ecoe"" (64) 

where c = eci/2. □ 

Let us denote the transformed variables by = (p^,q^)-^. The transformed variables are 
defined by ^'^(Ze) = z. Wc can then obtain an estimate for K(j)^{t)) in these coordinates by 



|if(p,(T))-i^(p,(0))| < |r| 



|:if(Pe(r)) 



< \r\\{K,H}\, 

< \r\er^^^\{K,R}\, 

< de-"^^^', (65) 

for \t\ < e'^l'^^ and some constant d > 0. This means that ||pj(t)|P stays almost constant for very 
long time intervals. 

Corollary 2.2. LeA us assume that the momenta p = q satisfy 

p(0) = -eJ2Vqy(q(0)) + Oie^) (66) 

at initial time t — 0, then 

p(r) = -eJ2Vqy(q(T)) + Oie^) (67) 

for all\T\< e"^^^. 

Proof. To first order in e the slow manifold is determined by p = — eJ2Vqy(q) and, in terms of 
the transformed coordinates, by p^ = C(e^)- Hence H66f) implies that K{p^{0)) = C(e^)- Finally 
(|65|l yields the desired result. □ 

This means that provided the system is within 0{e'^) of the geostrophically balanced state 
initially, it will stay there for exponentially long time intervals. It is not necessary for the system 
to be exactly geostrophically balanced for the the result to hold. 

If Pe(0) = 0, {i.e., we start on the slow manifold) then there are no rapid oscillations due to 
Pe, and we get the symplectic slow equation given in the following corollary: 

Corollary 2.3. //pe(0) = then 

J2q. = £VqG(0, q,) + 0{e~^^^'), (68) 

for aU\T\< 

Proof If Pe(0) = 0, then theorem [Q says that Pe{t) stays exponentially small for exponentially 
long time intervals and so 

q, = £VpG(pe,q,)|p.=o + 0(e-=/2^'), (69) 

for all |r| < e^/^"^. Wc can obtain the slow equations in symplectic form by examining the pg 
equation and ignoring exponentially small terms, i.e., 

= p, = J2VpG(pe,qe)|p.=0 - VqG(pe,qe)|p.=0, (70) 

and then substituting equation (|69|l to make 

J2qe = £VqG(p„qe)|p,=o = £VqG(0,q,), (71) 
up to terms exponentially small in e. □ 



We note that is in canonical form, i.e., the equations are Hamihonian with structure matrix 
Jisg = Jj and Hamiltonian 7Jisg(qe) = eG(0,qj). When the system is not initiahsed on the slow 
manifold then the dynamics is dominated by this slow equation but also contains fast components 
(which remain approximately of the same magnitude for exponentially long time intervals). 

In section 12.61 we will conduct a simple experiment that will allow us to verify the above 
estimates numerically. The particular set-up will be chosen such that p = Pe at the initial time 
r = and final time t = T . Then |_ftr(p(0))— X(p(T))| should go to zero exponentially fast as e ^ 
(discounting numerical round-off errors). We will show in scction r2.5l that such exponentially small 
terms can indeed be observed using a symplectic time stepping method (see, e.i?., |HLW02I ILR,05j 
for a discussion of symplectic time stepping methods for Hamiltonian ODEs). 

In the appendix we calculate the slow equation H()8|l to second-order by applying the iterative 
algorithm and obtain 

J2Cie = £VqF(q,) - YVq|lVqF(qe)|l2 + 0{e^). (72) 

Note that the significant terms in (|72() are equivalent to the 'large-scale semi-geostrophic' equation 
(j35|l . In other words, 1)68(1 provides a higher order generalisation of the 'large-scale semi-geostrophic' 
theory in the context of our simple model system. More precisely. Hamiltonian normal form theory 
naturally leads to A = —1/2. 

2.4 Extension to non-autonomous systems 

To recover results for equations ©-(O we need to extend the results from section to systems 
with time-dependent potentials where the Hamiltonian takes the form 

H„{^,T)^K{p)+eV{ci,T). (73) 

We apply the standard technique of extending the phase space by two extra conjugate variables s 
and e, and write a new Hamiltonian 

i/o(z, e, s, r) = i^(p) + el^(q, .s) - e, (74) 

so that the equations of motion for e and s are 

dV 

s = 1, e = e— (q, s), (75) 

and we recover the original equations. The new Hamiltonian Hq is now autonomous and in a suit- 
able form to apply theorem 12 . II with minor modifications (to account for the two extra variables), 
i.e., there exists a near-identity symplectic change of coordinates (p, q, e, s) — > (pe, q^, e^, s^) such 
that K{Tp^) is nearly preserved over exponentially long time intervals. 

We will make use of this extension to non-autonomous systems again in section [3. 21 

2.5 Numerical methods and backward error analysis 

The Hamiltonian equations H14() - (|15|l can be solved numerically using a symplectic method. As an 
example, we give the following splitting method. We rewrite the Hamiltonian (|23|l as a sum 

i?o - \k + eV+^K (76) 

and note that each entry has an exact flow map, which we denote by ^r,K/2 and ^r.eV, respectively. 
A second-order symplectic method is now given by the composition of these flow maps |LR,05j . i.e.. 



Z"+1 = 7\/Ar(z"), Mat = $Ar/2,A' O $Ar,eV O 'J'Ar/2,X- 



(77) 



Using backward error analysis, it can be shown that this method is equivalent to the exact 
solution of a modified Hamiltonian problem up to terms exponentially small in At |BG94I IRei99aj . 
More specifically, there exists a modified Hamiltonian 

i/Ar(z) =i?o(z)+£Ar2p(z,Ar) (78) 

such that 

||A/a.(z) - ^Ar.H^MW < ci e-^^/'^^ (79) 

where ci,C2 > arc appropriate constants. See |HLW02I ILROSj for a survey of backward error 
analysis results for Hamiltonian systems. 

The modified Hamiltonian H78(l allows us to discuss the preservation of geostrophic balance 
under the time-stepping method H77(l . All one has to do is to apply the normal form transformation 
from the previous section to the modified Hamiltonian (|78|) and to generalise theorem 12.11 In 
particular, the following theorem implies the near conservation of geostrophic balance under the 
numerical method 1)77(1 over exponentially long times. 

Theorem 2.4. Consider a symplectic integrator 

z"+i = MAt(z"), z" = ((q")^, (p")^)^, (80) 

such as j7?| j. Then, for a sufficiently small time step At, there exists a symplectic change of 
coordinates ^'e,Ar such that 

WKXr ° ^^Ar ° *e,Ar (Ze) - ^Ar,H^M)\\ < ^1 e"'^^/^^ (81) 

where = 'I'^Arl^) denotes the transformed variable, is the exact flow map of a Hamil- 

tonian system with transformed Hamiltonian 

Hat = Hat O *e,Ar =K + eG Ar + er'^'^/'RAr, (82) 

Hat is the exponentially accurate modified Hamiltonian for Mai fi-e-, j7i^ holds), and di > 0, 
i = 1,2,3, are appropriate constants. The transformation is chosen such that 

{K, Gat} = 0. (83) 

Furthermore, as Ar — s- 0, 

Gat G, Rat R, (84) 
where G and R are given in theorem \2.1\ 

Proof. The proof combines normal form estimates with backward error analysis for symplectic 
methods as developed in |R,ei99bj for adiabatic invariants in highly oscillatory mechanical systems. 
In particular, backward error analysis implies that the trajectory of the symplectic method Mat 
stays Ar-exponcntially close to the exact solution of a system with a modified Hamiltonian Hat 
consisting of iJo perturbed by an ©(eAr^) correction, p > 1 the order of the method. We then 
apply theorem 12. II to the modified Hamiltonian Hat- CH 

Thcorcm l2.4l imDlies that if p° is initialised at zero in the transformed coordinates, then p" will 
stay exponentially close to zero in the transformed coordinates over exponentially long times. This 
statement agrees with general results on the preservation of adiabatic invariants under symplectic 
time-stepping methods. See, for example, |Rei99b| and |LR05| . We will make use of this property 
of symplectic methods in the following section. 

Another immediate (more practical) conclusion from theorem 12.41 is the following 

Corollary 2.5. Numerical trajectories computed with a symplectic method such as jy?| ) satisfy 

p" = -£J2Vql^(q") + Oie^), (85) 

over exponentially many time-steps n provided the conditions of thcorem \2.4\ are satisfied and 1^85]) 
holds at n ^ 0. 



Proof. The proof is similar to that of corollary 12. 21 



□ 



2.6 A numerical experiment 

Wc consider the particular potential energy function 

^(q) ^<l- + \ cxp(-(g' + ql)), q = {qx,qyY, (86) 

and initial conditions 

g,(0) = 0, 9^(0) < -10, P:.(0) = 0, Py{Q) = e, (87) 

i.e., T^(q(0)) ~ 0. Simulations are run for long enough such that qy{T) > 10 at final time T. The 
nonlinear dynamics reduces to a linear system 

= 0, qy = e (88) 

at both T = and t ^ T, i.e., to a system in perfect balance. More precisely, according to the 
theoretical results of sections 12.31 & 12.51 this statement should be true at final time up to terms 
exponentially small in e. 



kinetic energy drift 
exponential fit 
total energy drift 




e 



Figure 1: Drift in total energy, kinetic energy, and exponential fit with 4 exp(— 0.92/e) as a function 
of e. 



We conduct a numerical experiment to verify this claim. We implement the symplectic method 
(|77|l with step size At = and vary e in the interval r G [1/15,1]. The drift in total energy 
AE = |i?(q(0),p(0)) - i?(q(r),p(T))| and in kinetic energy AK = \K{p{0)) - KipiT))\ is 
monitored as a function of e and the results are displayed in figure ^ The drift in total energy is 
solely caused by the symplectic time stepping method. The choice of At ensures that this drift is 
orders of magnitude smaller than the computed drift in kinetic energy. In other words, finite step 
size effects on the generation of unbalanced motion can be neglected. The anticipated exponential 
dependence of the drift in kinetic energy is nicely confirmed and 

AA" < 4 exp(-0.92/e) (89) 



provides an excellent boimd confirming theorem 12 .41 with ^3 = 0.92. 



3 Particle methods for the shallow-water equations 

Throughout scction|21 wc have assumed that the fluid depth /i is a given function of space and time 
and derived equations of motion for a single 'fluid parcel'. In this section, we consider a numerical 
approximation for fi of the form 

N 

/^(t,x) -^m,7A(||x-q,(t)||) (90) 

i=l 

in terms of N moving 'fluid parcels', where q_i{t) G M.'^ denotes the location of the ith fluid parcel 
at time t with mass rrii and shape function ')/'('') ^ 0. The approximation (|90|l provides the typical 
starting point for a particle method such as smoothed particle hydrodynamics (SPH), which was 
first proposed in |Luc77l [GMTfj for general fluid dynamics and for the SWEs CJ-® in |Sal83j . 
Each 'fluid parcel' moves under the Newtonian equations of motion 

d B 
^°di^' -^zPi - ^VxM(t,x)|x=q., (91) 

= P, (92) 

The equations ()9()(I - H92|I form a closed set of equations and provide an approximation to ©-(CI)- 
(See jFG!R,02l IFR.031 IFR,04I IBHR.05j for a numerically robust implementation of a particle method 
for the SWEs and their geometric properties.) 

After setting B = Ro = e and a rescaling of time, equations (|91|) - (|92|l become 

Pi = J2Pi - eVxM(x,i)|x=qi, (93) 
= Pi- (94) 

Without restriction of generality, we may also assume that all 'fluid parcels' carry a constant mass 
rrii = S. The equations (|9l)|l . (|93|I - H94|I are Hamiltonian with Hamiltonian function 

N . N 



=1 



structure matrix 



phase space variable z = (q^,P"^)"^ e R^^, and q = (qf , . . . ,q|^)^, p = (p^, . . . ,p^)"^. Note 
that a finite fluid depth fi implies that the particle masses rrii = S approach zero as the number of 
particles N — > oo. More precisely, S N const, as N oo. 



3.1 Preservation of geostrophic motion 

We first note that the Hamiltonian system defined by H95|) - (|96|) with = 1 fits exactly into 
the framework of section 2. Furthermore, as first observed in |Cot04| . the normal form theory 
developed in section 2 can be generalised to > 1 without much difficulty. The key idea goes 
back to an observation in |BC!G87j . Specifically, it turns out that the normal form theory for a 
system with many degrees of motion with identical fast frequency proceeds much along the lines 
of the theory for a single fast degree of motion. Most importantly, the exponential estimates do 
not depend on the number of degrees of freedom and we obtain: 

Theorem 3.1. We consider a fixed number of particles N . Let the shape function ip be analytic in 
some compact subset of C?^ containing the solution trajectory. Then there exists a near-identity 
change of coordinates such that 



H = i/sph o^,_= K + eG + e~''^ R 



(97) 



where 

{G, K} = Q (98) 

and c > is some constant. 

Proof. We apply theorem 12. ll in slightly modified form taking into account that > 1. We also 
note that a bound on ijj immediately implies a boimd on the potential energy V via \V{c{)\ < 
SN\^P\oo. □ 

The estimate (|65f) hence also applies to the particle approximation 193|) - (|94|l and geostrophic 
balance is maintained over exponentially long time intervals. Furthermore, theorem l2.4l and corol- 
lary!^] immediately carry over to the particle method. Note that H85|l gets replaced by 

= -£ J2Vxm(x, t)|x=q. + 0{e') (99) 

for i — 1, . . . ,N. Recall that we consider the limit e ^ for N fixed. 

It is important to keep in mind that theorem 13.11 does not apply to most practical implemen- 
tations of particle method since the basis functions are typically not analytic. In those case the 
exponentially small term in (|97|l needs to be replaced by a polynomial expression of the form (ce)*^, 
where the integer k > 2 and the constant c > depend on the smoothness of the basis function ip. 
This is the case, for example, for the Hamiltonian-Particle-Mesh (HPM) method (see |FGR,02| 1. 
which will be used in the following section to conduct a number of numerical experiments. 



3.2 Numerical experiments 

In this section, two numerical examples are given which illustrates theorems l2.4l and l3.1l for practical 
implementations of a particle method. 



3.2.1 Exchange of kinetic energy 

To show that particles may exchange energy as long as the total kinetic energy is preserved, 
consider a model consisting of 2 interacting particles moving on a predefined background height 
field 7i(T, x). In this experiment, we use the Hamiltonian Particle-Mesh (HPM) method for spatial 
discretisation (see |FGR02| ') with the total fluid depth at a grid point Xm„ given by 

2 

Mmn(T) = y^V'mn(qi(T)) + Jl{T,X.rnn) ■ (100) 

1=1 

The basis functions ipmn are normalized cubic B-splines centered about grid points Xm„. 

We multiplied the potential energy by an analytic function (/(r) with g{T) ^ as t ^ ±00 so 
that the Hamiltonian takes the form 

^ 1 e ^ / 2 \ 

Hi,p„,{z) ^^-\\p^\\'^ + -g{T)^^t(j^n{cii) ^ V„m(qj) + 27l(r,x™„) . (101) 

i—1 mn i—1 V.?"-'^ / 

The reason for the introduction of the function ^(t) is that when g{T) « at the beginning and the 
end of the simulation, the Hamiltonian is already in normal form and we can measure the amount 
of fast rotational energy exactly (this technique was first introduced in |BF99j in an experiment to 
measure the exchange of oscillatory energy between two colliding molecules). 

A graph showing kinetic energy against time during the experiment, and a plot of the trajec- 
tories of the two particles is shown in figure |2| This example illustrates that, although the total 
kinetic energy K is almost invariant, the particles exchange energy through the interaction poten- 
tial as shown in figure 01 This exchange is permitted because the frequencies of oscillation of the 
two particles are the same for e = 0. 
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Figure 2: Plots showing particle trajectories for a system of 2 interacting HPM particles moving 
on a prescribed background height field. The initial positions are marked with a and the end 
positions are marked with a '+'. The kinetic energy for each particle is shown in figure |31 
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Figure 3: Plots showing kinetic energies for 2 interacting HPM particles moving on a prescribed 
background height field. The total kinetic energy returns at the end of the experiment to a value 
very close to the starting value. However the 2 particles do exchange kinetic energy. 



3.2.2 Preservation of geostrophic balance under the HPM method 

As already eluded to, theorem 13 . II docs not directly apply to the HPM method |F(Tri,02j since the 
basis functions are not analytic. However, one would nevertheless expect 'good' preservation of 
geostrophic balance in the semi-geostrophic scaling limit. To test this hypothesis numerically, we 
repeated the shear flow instability simulation of |FG!R02j under slightly modified initial conditions 
such that the associated Rossby and Burger number satisfy Ro Bu sa 0.1. The HPM discretized 
shallow water equations are simulated using TV = 262144 particles over a domain (a;, y) G [0, 27r]^, 
a smoothing length of a = 0.2015, and a time step of At = 1/36. Note that one unit of time 
corresponds to one day and that the Rossby radius of deformation corresponds to « 0.5 in 
dimensionlcss variables. The initial momenta are in geostrophic balance, i.e., 

p(0) = -eJ2ArVqy(q). (102) 
fluid height at initial time fiuid heigiit at finai time 
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Figure 4: Drift in total energy, kinetic energy, and exponential fit with 4 exp(— 0.92/e) as a fimction 
of e. 

Simulation results can be found in figure 21 where we display the initial and final layer depth 
{t = 15 days). We also monitor the total kinetic energy K{p), the total kinetic energy in the 
ageostrophic momentum components 

Pag = P - Pgs, Pgs = -eJ2wVqF(q), (103) 

i.e., K{piig), and the relative error in total energy i/hpm- We find that the relative error in energy 
remains very small. This is due to the symplectic nature of the time stepping procedure. We also 
observe that if(pag) remains nearly level and much smaller than the kinetic energy K{p). This 
indicates that geostrophic balance is well preserved throughout the simulation even though e sa 0.1 
is not very small and the basis functions 4'mn are not analytic. 

4 Summary and outlook 

We introduced the model system of a single fluid particle moving in a rotating shallow-water system 
and gave an exponentially-accurate normal form theorem valid in the semi-geostrophic scaling 



limit. This normal form theorem gives a coordinate change for which the kinetic energy stays 
almost invariant for very long time intervals. When the kinetic energy in transformed variables is 
initially zero (or close to zero) we showed that it is possible to derive an exponentially-accurate slow 
equation which is equivalent to the symplectic slow equation obtained by variational asymptotics 
in the semi-geostrophic limit. 

We showed that this result may be extended to numerical schemes for solving the SWEs based 
on particle methods by making use of backward error analysis, and then illustrated this with a 
numerical experiment designed to "capture" the exact amount of energy exchanged between the 
slow and fast dynamics. 

The extension of this work to the full rotating shallow- water PDE requires bounds 

to be obtained on the gradients of the height field. As the SWEs are known to develop shocks 
it seems unlikely that an exponentially-accurate result will be valid for long times unless it is 
possible to obtain some extra regularity when the fast dynamics has very little energy. It may be 
possible to obtain estimates for regularised equations where a smoothing operator is applied to the 
height gradient in the momentum equation, or for the family of a-regularised equations such as 
shallow- water-a |Hol99j . 

Finally we suggest one possibly fruitful application of this work in the calculation of the dis- 
persion of passive tracers using Lagrangian particles. Typically these calculations arc performed 
by interpolating a given (gridded) velocity field Umn{t) from the grid to the whole domain, and 
solving the system of ODEs 

^ =Y,4'mni(l^)u^nit), l=l,...,N. (104) 

mn 

This method can often experience problems with artificial clumping of particles when the trajec- 
tories are integrated over long time intervals. An alternative approach would be to solve equations 
with a given fluid depth (or geopotcntial) ji. The results of this paper show that if the 
system is integrated with a suitable symplectic method then the tracer particles will keep balanced 
trajectories in the semi-geostrophic limit. 
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Appendix. Calculation of slow equation in symplectic form 
to second-order 

In this appendix we explicitly calculate the change of coordinates and the slow equation to 0{e^) 
using the iterative procedure given in section [2.31 

The Hamiltonian vector field given by the kinetic energy K = ||p|p/2 is 

p = J2P, q = p, 

and the associated time-r flow map ^t,k is 

p e'^^'^p, 

q ^ J2(l-e^^")p + q. 

First we need to calculate the transformed Hamiltonian after the first-order transformation. 

= Ho + e{Ho, Fi} + ^{{i^o, i^i}, F^} + 0{e^), 

= K + eV + s{K, Fi} + s^V, 51} + ^{{A', FJ, Fi} + ©(e^), 
- - 

= K + eV + —{V + V,Fi} + Oie^), 
where we have used {K,gi} ~ V — V, {K, V} = 0, and 

" ^^^^(•^2(l-e^^^)p + q), 

as well as 

Fi(z) ^ ^ I drr(y-F)o$,_^(z), 
^ Jo 

- drr I y(j2(l-e'^-^)p + q) dr' F ( J2 (l - e'^^(^+^') 



p + q 



Recall that z = (p^, q^)"^ and T = 27r. 

The first-order slow equation is given by: 

J2q = eVq^ (p, q) |p=o = eVqy(q). 

To calculate the slow equation at the next order, we write the transformed Hamiltonian after the 
first-order transformation as 

H2 = K + eV + e^Ri+0{e^), 



where 



^1 = 2^V + V,F,}. 



Note that for a general function /(p,q), 

Vq/(P, q)|p=o = ^ ^ dr Vq {/ (e'^^^p, J2 (l - e^^^) p + q) } |p=o = Vq/(0, q), 

so all that remains is to calculate Ri and the gradient with respect to q for p = 0. 
We have 

{V, Fi} - VpV ■ JzVpFi + Vql7 • VpFi - VpV ■ VqFi, 

and 

= VqF- Vpi^i. 

Now we need to calculate these various gradients: 

Vqi^i(z)|p=o = ^/ rfTTVq(y- V)o$,,^(z)|p=o, 



= Vq(V^(q) - V^(q)) = 0, 



and 



VpFi(z)|p=o 



J2 (1 - e^^^) VqF ( J2 (1 - e'^^") p + q) |p=o - 
dr' J, (1 - e;^^('^-+'^-')) Vqy (J2 (1 - e'^^^) p + q) |p=o 
drr .h{l-e''^^)- £ dr'J,[l-e'^^-^^^+^^'^) ^ Vq^q), 



dTT-^e^^^VqV^(q), 



-Vqy(q). 



dr 



We also obtain 



Vqy(z)|p=o = ^/ cfTVqy(J2(l-e^^-)p + q))|p=o, 



Vqy(q). 



and 



Vpy(z)|p=o = ^/ J2VqF(J2(l-e-'^^)p + q) Ip^o, 



Hence we get 
and 



= -J2Vqy(q). 

{V, i^l}|p = - - (J2Vql/) • ( J2Vqy) - Vq^ • Vql/ = 0, 
{V,Fl}\p = = -Vq^- Vq^ = -||Vqy|P, 



and so i?i = ^^llVqFjp. The slow equation in transformed coordinates is then 

We note that this is the same slow canonical equation obtained using variational asymptotics in 
|(T()05al IGO05bj and. in particular, yields the 'large-scale semi-geostrophic' equation l|35(l . 
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